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Abstract 

In a recent preprint ( arXiv: 121 1.4285^ 1) [T4] we addressed the 
problem of constructing reduced models for time-dependent systems 
described by differential equations which involve uncertain parame- 
ters. In the current work, we focus on the construction of reduced 
models for systems of differential equations when the initial condition 
is uncertain. While for both cases the reduced models are constructed 
through the Mori-Zwanzig formalism, the necessary estimation of the 
memory parameters is quite different. For the case of uncertain initial 
conditions we present an algorithm which allows to estimate on the 
fly the parameters appearing in the reduced model. The first part of 
the algorithm evolves the full system until the estimation of the pa- 
rameters for the reduced model has converged. At the time instant 
that this happens, the algorithm switches to the evolution of only the 
reduced model with the estimated parameter values from the first part 
of the algorithm. The viscous Burgers equation with uncertain initial 
condition is used to illustrate the construction. 



1 Introduction 



The problem of quantifying the uncertainty of the solution of systems of par- 
tial or ordinary differential equations has become in recent years a rather 
active area of research. The realization that more often than not, for prob- 
lems of practical interest, one is not able to determine the parameters, initial 
conditions, boundary conditions etc. to within high enough accuracy, has 
led to a flourishing literature of methods for quantifying the impact that this 
uncertainty imposes on the solution of the problems under investigation (see 
e.g. [Si mi [l2l [131 [m [E] ) . However, despite the increase in computational 
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power and the development of various techniques for uncertainty quantifica- 
tion there is still a wealth of problems where reliable uncertainty quantifica- 
tion is beyond reach. The main reason behind the inadequacy is the often 
high dimensionality (in probability space) of the uncertainty sources. When 
this uncertainty is coupled with the fact that for practical problems, even 
solving the corresponding equations for one value of the uncertain parameter 
(initial condition, boundary condition, . . .) can be very expensive, it results 
in the uncertainty quantification problem being a rather formidable task. 
One way to address this problem is to look for reduced models for a subset 
of the variables needed for a complete description of the uncertainty. 

In the current work, we are concerned with the construction of reduced 
models for systems of differential equations that arise from polynomial chaos 
expansions of solutions of a PDE or ODE system. In [14] we focus on 
the construction of reduced models when there exists uncertainty in the 
parameters of the original system. In the current work, we focus on the 
case that the given PDE or ODE system has uncertain initial condition. 
Our goal is to construct a reduced model for the evolution of a subset of 
the polynomial chaos expansions that are needed for a complete description 
of the uncertainty caused by the uncertainty in the initial condition. There 
are different methods to construct reduced models for PDE or ODE systems 
(see e.g. 011] and references therein). We choose to use the Mori-Zwanzig 
(MZ) formalism in order to construct the reduced model [2l[3]. 

The main issue with all model reduction approaches is the computation 
of the memory caused by the process of eliminating variables from the given 
system (referred to as the full system from this point on) [3]. The memory 
terms are, in general, integral terms which account for the history of the 
variables that are not resolved. In principle the integrands appearing in 
the memory terms can be computed through the solution of the orthogonal 
dynamics equation [3]. However, the solution of this equation is usually 
very expensive. As we did in [14], here too, we utilize a Markovian refor- 
mulation of the MZ formalism which allows the calculation of the memory 
terms through the solution of ordinary differential equations instead of the 
computation of convolution integrals as they appear in the original formu- 
lation. For the example presented in [11], the parameters appearing in the 
reduced model could be estimated without having to solve the full system. 
In the current work, this is not possible. However, we present an algorithm 
which allows the estimation of the necessary parameters on the fly. This 
means that one starts evolving the full system and use it to estimate the 
reduced model parameters. Once this is achieved, the simulation continues 
by evolving only the reduced model with the necessary parameters set equal 
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to their estimated values from the first part of the algorithm. 

Section [2] presents a brief introduction to the MZ formalism for the con- 
struction of reduced models of systems of ODEs. In Section [3] we develop 
the Markovian reformulation of the MZ formalism and show how one can 
estimate adaptively the parameters appearing in the reduced model. Sec- 
tion [3] applies the reformulation of MZ presented in Section [3] to the viscous 
Burgers equation with uncertain initial condition. Finally, in Section [5] we 
discuss certain directions for future work. 



2 Mori-Zwanzig formalism 

We begin with a brief presentation of the Mori-Zwanzig formalism [21 [3]. 
Suppose we are given the system 

^ = R{t,um (1) 

where u = {{uk}), k G H Li G with initial condition u(0) = uq. Our goal is 
to construct a reduced model for the modes in the subset H. The system 
of ordinary differential equations we are given can be transformed into a 
system of linear partial differential equations 

^"^^ L<Pk, Muo,0) = uok,keHuG (2) 



dt 



where L = ^keHuG ^i('^o)d^i ■ "^^^ solution of ^ is given by Uk{uo,t) = 
(pk (uq J t) ■ Using semigroup notation we can rewrite ([2]) as 

d 
dt 

Suppose that the vector of initial conditions can be divided as uq = (no, no), 
where no is the vector of the resolved variables (those in H) and no is the 
vector of the unresolved variables (those in G). Let P be an orthogonal 
projection on the space of functions of no and Q = I — P. 
Equation ([2]) can be rewritten as 



d_ 

dt 



e^^uok = e^^PLuok + e^^^QLuok + [ e^^''^^ PLe'^^QLu^kds, k e H, 

Jo 

(3) 

where we have used Dyson's formula 

gtL ^ ^tQL ^ /■* e{t-^)ip2,e"^-^ds. (4) 
Jo 
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Equation Q is the Mori-Zwanzig identity. Note that this relation is exact 
and is an alternative way of writing the original PDE. It is the starting 
point of our approximations. Of course, we have one such equation for each 
of the resolved variables Uk,k G H. The first term in ([3]) is usually called 
Markovian since it depends only on the values of the variables at the current 
instant, the second is called "noise" and the third "memory". 
If we write 

e^^^QLuok = Wk, 
Wk{uo,t) satisfies the equation 

§lWk{uQ,t) = QLwk{uo, t) 
Wkiuo,0) = QLxk = Rkiuo) - {PRk){uo). 

If we project ([5]) we get 

d 

P—Wk{uo,t) = PQLwk{uo,t) = 0, 
since PQ = 0. Also for the initial condition 

Pwk{uo,0) = PQLuok = 

by the same argument. Thus, the solution of ([5]) is at all times orthogonal 
to the range of P. We call ([5]) the orthogonal dynamics equation. Since 
the solutions of the orthogonal dynamics equation remain orthogonal to the 
range of P, we can project the Mori-Zwanzig equation ([3]) and find 

^Pe'^uok = Pe'^PLuok + P [ e^'~'^^ PLe'^^QLuokds. (6) 
ot Jo 

3 Finite memory 

In this section we describe a reformulation of the problem of computing the 
memory term which does not use the orthogonal dynamics equation. We fo- 
cus on the case when the memory has a finite extent only. The case of infinite 
memory is simpler and is a special case of the formulation presented below. 
Also, the current reformulation allows us to comment on what happens in 
the case when the memory is very short. 

Let wok{t) = P Jq e^^-'^^PLe"^^QLuokds = P e'^PLe^^-'^^^QLuQkds, 
by the change of variables t' = t — s. Note, that w^k depends both on t and 
the resolved part of the initial conditions uq. We have suppressed the uq 
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dependence for simplicity of notation. If the memory extends only for to 
units in the past (with to <t,) then 



I t-to 

The evolution of wok is given by 

dwok 
dt 

where 



{t) = P [ e'^PLe^^-'^^^QLuokds. 

Jt-tn 



= Pe^^PLQLuok - Pe(*-*°)^PLe*o^^QLnofe + wik(t), (7) 



wik{t) = P [ e'^PLe^^-'^^^QLQLuQkds. 
Jt-to 



-to 

To allow for more flexibility, let us assume that the integrand in the formula 
for wik{t) contributes only for ti units with ti < to- Then 



{t) = P f e'^PLe^^-'^^^QLQLuokds. 
Jt-t^ 



Wlk 

it-ti 

We can proceed and write an equation for the evolution of wik(t) which 
reads 



dt 

where 



Pe^^PLQLQLuok - Pe^^-^^^^PLe^^^^QLQLuok + W2k{t), (8) 



W2k{t) = P [ e'^PLe^^-'^'^^QLQLQLuQkds. 
Jt-ti 

Similarly, if this integral extends only for t2 units in the past with t2 <ti, 
then ^ 

W2kit) = P [ e'^PLe^^-'^^^QLQLQLuokds. 

Jt-t2 

This hierarchy of equations continues indefinitely. Also, we can assume for 
more flexibility that at every level of the hierarchy we allow the interval of 

integration for the integral term to extend to fewer or the same units of 
time than the integral in the previous level. If we keep, say, n terms in this 
hierarchy, the equation for wi^n-i)k{t) will read 

^^^^^^ = Pe'^^PLiQir-^QLuok- (9) 

Pg(t-t„_l)Lp^gt„_lQL(g^)n-lg^^^^ ^ ^^^^^^ 
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where 

Wnk 



Note that the last term in ([9]) involves the unknown evolution operator for 
the orthogonal dynamics equation. This situation is the well-known closure 
problem. We can stop the hierarchy at the nth term by assuming that 

Wnk{t) = 0. 

In addition to the closure problem, the unknown evolution operator for 
the orthogonal dynamics equation appears in the equations for the evolu- 
tion of WQk{t), • • • , W{n-i)k{i) through the terms Pe^*~*°^^PLe*°^^QLuok, ■ ■ ■ 
Pe(*-*o)^PLe*o'3^(QL)"-iQLuofc respectively. 

We describe now a way to express these terms involving the unknown 
orthogonal dynamics operator through known quantities so that we obtain 
a closed system for the evolution of wok{t), . . . , W(^n-i)ki't)- 

Since we want to treat the case where to is not necessarily small, we 
divide the interval [t — tQ,t] in uq subintervals. Define 



(1) 



(t) = P f e'^PLe'^'-'^'^^QLuokds 

Jt~Atn 



Ato 
t-Ato 



^Okit) = P I e'^PLe^'-'^^^QLuokds 



Jt~to 



where no Ato = and WQk{t) = J2i^li''^ok(^)- Similarly, we can define the 
quantities w^j} (t), . . . , w^^^\t) 

u>f^)(t) = P / e'^PLe^^-'^'^^QLQLuokds 
Jt-Ati 
rt-Ah 

wffjit) = P / e'^PLe^'-'^^^QLQLuokds 



t-{ni"l)Ati 
t-t, 



w^;i'\t) = P / e'^^PLe^'-'^'^'-QLQLuQkds, 
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where niAti = ti and wik{t) = Y17=i ''^il ■ ^ similar fashion we can 
define corresponding quantities for aU the memory terms up to = 

In order to proceed we need to make an approximation for the integrals 
over the subintervals. 

3.1 Trapezoidal rule approximation 

We have 



Jt-Ato 
from which we find 

and from ([7]) 

(2) 

Similarly, for w^^ (t) we find 

(2) 



dt V 



At. /^ofe - ^Pe'^'PLQLuok + w\^^{t) + 0{{Ato) 



Co 

In general, 



-j-i 

+ ' 



+ <^^(t) + 0((Ato)') for i = l,...,no. 

(10) 
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Similarly, 



dw 



dt \Ati 

+ 



+ w^lit) + 0((Ati)2) for i = l,...,ni 



-j=i 
dw^"^ 



dT' - -yAt^^j'^ln^Dk^t) + {-iy+'2Pe''^PL{QLr-'QLuok 



+ 



+ 0((Ai„_i)2) for i = l,...,n„_i. 



fll) 



By dropping the 0((Ato)^), • • • ,0((At„_i)^) terms we obtain a system of 
riQ + ni + . . . + Un-i differential equations for the evolution of the quan- 
tities w'^Q^{t), . . . This system allows us to determine the memory 

term tt;ofc(i) = P Jq e^^'"^^^ PLe'^'^^QLuokds. Since the approximation we 
have used for the integral leads to an error 0(Ai)2, the ODE solver should 
also be 0(At)^. We have used the modified Euler method to solve numeri- 
cally the equations for the reduced model. 

Note that the implementation of the above scheme requires the knowl- 
edge of the expressions for Pe^^PLQLuok-, ■ ■ ■ , Pe*'" PL{QL)"-~^QLuQf^. Since 
the computation of these expressions for large n can be rather involved for 
nonlinear systems (see Section H]), we expect that the above scheme will 
be used with a small to moderate value of n. Finally, we mention that the 
above construction can be carried out for integration rules of higher order 
e.g. Simpson's rule. 



3.2 Estimation of the memory length 

The construction presented above relies on an accurate determination of the 
memory lengths Ato, Ati, . . . .At„_i. We present in this section a way to 
estimate these quantities on the fly. This means that we start evolving the 
full system, use it to estimate Ato, Ati, . . . .Atn-i and then switch to the 
reduced model with the estimated values for Ato, Ati, . . . .At„_i. 

For simplicity of presentation we assume that we evolve only WQkit) and 
use only one subinterval to discretize the time integrals, i.e. AtQ = to- The 
reduced model reads 
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^ = Pe''^PLuok + wok{t) (12) 

^ = 2Pe*^PLQLnofc - r^ofe(i) (13) 
at iQ 

for k € F. We can solve (fT3]) formally and substitute in to get 

^ = Pe'^PLuok + t e'^'^^'-'hPe'^PLQL 
at Jo 



uokds (14) 

where Aq = 2/to- Recall that, for the resolved variables, we have from the 
full system 

Pe^^PLuok + Pe^^QLuok. (15) 



dt 

We would like to estimate the memory decay parameter to so that the re- 
duced equation ([H]) for Uk reproduces the behavior of Uk as predicted by 
the full system (jlSp . We can do that by requiring that the evolution of some 
integral quantity of the solution is the same when predicted by the reduced 
and full systems. 

We begin by discretizing the integral term in (114p . Suppose that we are 
evolving the full system with a step size 5t, where t = ritSt (note that nt 
increases as t increases). If we discretize the integral with the trapezoidal 
rule we find 

^ = Pe'^^PLuok (16) 

nt-l ^ 

+ [A(t,no) + 2 5] e-^"(*-^'^*)A(i5t,^io) + e-^"Vfc(0,no)]y 
i=i 

where fk{jSt,uo) = 2Pe^^^^ PLQLuQk for j = 0, ...,nt. The quantities 
fk{jSt,uo) can be computed from the full system. 

There is freedom in the choice of the integral quantity whose evolution 
the reduced model should be able to reproduce. For example, we can use 
Z]fcG-F squared I2 norm of the resolved variables. If we use this 

integral quantity, then from ()16p and (jlSp we find that the unknown param- 
eter to must satisfy 

^ 2Re{Ik{t, to)ul{t)} = 2i?e{Pe*^QLnofc4(t)}, (17) 

k£F keF 
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where 



Ik{t,to) = [A(t,no) + 2 ^ e-^°(*-^'5*)/fc(j5t,no) + e-^»*/fc(0,no)]- 

and Re{-} denotes the real part. 
Let y = exp[— Ao5t]. Then, 

nt-l e, 

Ik{t,to) = [A(t,no) + 2 ^ 2/"*-Vfc(i5i,^o) + y"7fc(0,tio)]y. (18) 

j=i 

With this identification, equation (jl7p becomes a polynomial equation for 
y with ?/ € [0, 1]. It is not difficult to solve equation (fTTl) with an iterative 
method, for example Newton's method. For the numerical results we present 
in SectionlU Newton's method converged to double precision accuracy within 
4-5 iterations. After an estimate y has been obtained, we can find the 
estimate to of (recall Aq = 2 /to) from 

2St , , 

to = -r-r. 19 
Iny 



3.2.1 Determination of optimal estimate 

For each time instant t we can obtain through equations (jl7p and (llQh . an 
estimate t^it) for to- Thus, the most important issue that we have to address 
is that of deciding which is the best estimate of to. In other words, at what 
time tf should we stop estimating the value of to so that we can use the 
estimated value to(tj) to evolve the reduced model from then on. 

We define e(t) = max |y'(t + bt) — y'(t)|. The quantity e(t) monitors 

/e[l,nt] 

the convergence of not only the value of the estimate y as a function of the 
time t, but of the whole function e--^o(i-s)_ ideally, e(t) converges to zero 
with increasing t. That will be the case if the approximation of the memory 
term only through Pe^^ PLQLuokr is enough (see ([12]) -([13])). However, this 
will not always be the case. If keeping Pe^^PLQLuQkr is not enough, then 
e(t) will decrease with increasing t up to some time tmin when it will reach 
a nonzero minimum. After that time, it starts increasing. This signals 
that keeping only Pe^^ PLQLuQkr is not enough to describe accurately the 
memory. 

In order to proceed we have two options: (i) construct a higher order 
model and (ii) identify tf = tmin and thus to(tj) = toitmin)- Results for 
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higher order models will be presented elsewhere (see also discussion in Sec- 
tion [5]). In the numerical experiments we present in the next section we have 
chosen to(i/) = io(imm)- Note that the procedure just outlined allows the 
automation of the algorithm. This means that there is no adjustable reduced 
model parameter that needs to be specified at the onset of the algorithm. 

We are now in a position to state the adaptive Mori-Zwanzig algorithm 
which constructs a reduced model with the necessary memory term param- 
eter to estimated on the fly. 

Adaptive Mori-Zwanzig Algorithm 

1. Evolve the full system and compute, at every step, the estimate to(^)- 
Use estimates of to from successive steps to calculate e(t) = max 

«G[l,nt] 

St)-y\t)\. 

2. When e{t) reaches a minimum (possibly non zero) value at some in- 
stant tmin, pick to^tmin) as the final estimate of Iq. 

3. For the remaining simulation time, switch from the full system to 
the reduced model. The reduced model is evolved with the necessary 
parameter to set to its estimated value to(*mm)- 

This procedure can be extended to the computation of optimal estimates 
for ti,t2,---, i.e. when we evolve, in addition to WQkit), the quantities 
wik (t) , W2k {i), ■ ■ ■ ■ Results for such higher order models will be presented 
elsewhere. 



4 Burgers equation with uncertain initial condi- 
tion 

In this section we show how the above MZ formulation can be used for 
uncertainty quantification for the one-dimensional Burgers equation with 
uncertain initial condition. The equation is given by 

ut + uu^ = vuxxi (20) 

where v > Q. Equation ()20p should be supplemented with an initial condition 
u(x, 0) = uq{x) and boundary conditions. We solve (|20p in the interval [0, 2tt] 
with periodic boundary conditions. This allows us to expand the solution 
in Fourier series 

UN{x,t) = ^nfc(t)e*''^. 
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where F = [~-2^ '2 — "^^^ equation of motion for the Fourier mode 
becomes 

— = -— 2^ UpUq - Vk Uk- (21) 
p+q=k 

We assume that the initial condition ^0(2;) is uncertain (random) and can 
be expanded as uq{x,^) = (oq + ai^)fo(x) where ^ is uniformly distributed 
in [—1, 1] and vo{x) a given function. In the numerical experiments we have 
taken oq = ai = 1 and vq{x) = sinx. 

To proceed we expand the solution Uk{t,(,) for A; € F in a polynomial 
chaos expansion using Legendre polynomials which are orthogonal in the 
interval [—1,1]. In particular, we have that 

1' L,(0L,(e)de = ^5ii, 

where ivj(^) is the Legendre polynomial of order i. For each wavenumber k 
we expand the solution Uk{t,^) of in Legendre polynomials and keep 
the first M polynomials 

M-l 

Uk{t,0 ~ ^k^{t)U{i), where i ~ U[-l, 1]. (22) 

i=0 

Similarly, the initial condition can be written as uq{x,^) = sin x Ylii=o '^i^iiO 
since Lo{^) = 1 and Li(^) = ^. Substitution of (i22]l in (i24|l . use of the ex- 
pansion of the viscosity coefficient and of the orthogonality properties of the 
Legendre polynomials gives 

duk (t) ik 

= ^ ^ ^ ^ ^ ^ ^ fJ'pl'U'qmClmr ^ Ukr (^3) 
1=0 m=0 p-\-q=k 
p,q&F 

lor k e F and r = 0, . . . , M - 1. Also 

where the expectation E[-] is taken with respect to the uniform density 
on [—1, 1]. The expectation on the denominator of the expression for q^r 
is E[L'^{^)] = J^i L'^iO^d^ = 2fq:T) while the expectation on the numer- 
ator can be computed accurately using Gaussian quadrature with Legen- 
dre nodes. The Legendre polynomial triple product integral defines a ten- 
sor which has the following sparsity pattern: E[Li{(,)L„i{(,)Lr{(,)] = 0, if 
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l + m<r or l + r<m or m + r<l or l + m + r = odd [8j . Due to this 
sparsity pattern, for a given value of M only about 1/4 of the tensor 
entries are different from zero. 



4.1 MZ reduced model 



To conform with the Mori-Zwanzig formalism we set 



M-lM-l 



1=0 m=0 p+q=k 
p,qeF 



where u = {ukr} for k ^ F and r = 0, . . . , M — 1. Thus, we have 



for k ^ F and r = 0, . . . ,M — 1. We proceed by dividing the variables in 
resolved and unresolved. In particular, we consider as resolved the variables 
u = {ufcr} for k (z F and r = 0, . . . , A — 1, where A < M. Similarly, the 
unresolved variables are u = {u^r} foi" k ^ F and r = A, ... ,M — 1. In 
the notation of Section [2] we have H = i*" U (0, . . . , A — 1) and G = F \J 
(A, . . . , M — 1). In other words, we resolve, for all the Fourier modes, only 
the first A of the Legendre expansion coefficients and we shall construct a 
reduced model for them. 

The system (pi]) is supplemented by the initial condition uq = (no,uo)- 
We focus on initial conditions where the unresolved Fourier modes are set 
to zero, i.e. uq = (no,0). We also define L by 



To construct a MZ reduced model we need to define a projection operator P. 
For a function h{uo) of all the variables, the projection operator we will use 
is defined by P{h{u)) = P{h{uQ,UQ)) = /i(uo,0), i.e. it replaces the value 
of the unresolved variables uq in any function h{uQ) by zero. Note that this 
choice of projection is consistent with the initial conditions we have chosen. 
Also, we define the Markovian term 




M-l 




A-1 A-1 



1=0 m=Op+q=k 
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The Markovian term has the same functional form as the RHS of the fuU 
system but is restricted to a sum over only the first A Legendre expansion 
coefficients for each Fourier mode. 

For the the term PLQLuQkr we find 



PLQLuokr = 2 X 



M-l A-1 

mr 

l=A m=Op+q=k 
p,qeF 



(25) 



Finally, to implement any method to solve equation (jl7p for the estima- 
tion of to we need to specify the RHS of the equation (fT7|l . This requires 
the evaluation of the expression Pe^^QLuQkr- For the case of the viscous 
Burgers equation, we find 

A/-1 A-1 

Pe^^QLuQkr = 2( -) X] X] UplUgmClmr (26) 

/=A m=Op+q=k 
p,q£F 

., M-lAf-1 

tk V ■\ V V 

^ / ^ / ^ / ^ UplUq-mClmr • 
l=A m=Kp+q=k 

Note that since we restrict attention to initial conditions for which the un- 
resolved variables are zero and the projection sets the unresolved variables 
to zero, the quantity Pe^^QLuQkr can be computed through the evolution 
of the fuh system (i24|l . 



4.2 Numerical results 

In this section we present numerical results for the reduced model of the 
viscous Burgers equation with viscosity coefficient u = 0.03. The solution of 
the full system was computed with A'' = 196 Fourier modes {F = [—98, 97]) 
and the first 7 Legendre polynomials (M = 7). The first 7 Legendre poly- 
nomials were enough to obtain converged statistics for the full system. The 
full system was solved with the modified Euler method with 5t = 0.001. 

The reduced model uses = 196 Fourier modes but only the first two 
Legendre polynomials, so A = 2. It was solved using the modified Euler 
method with 5t = 0.001. The parameter to needed for the evolution of the 
memory term was found to be 0.3783 through the procedure described in 
Section 13.2.11 
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Time 

Figure 1: Evolution of the mean of the energy of the solution using only the 

first ^'^^^'^ J ^pcff^nrlrf^ "nnK^nnmialt; 




Time 

Figure 2: Evolution of the standard deviation of the energy of the solution 
using only the first two Legendre polynomials. 
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Figure [T] shows the evolution of the mean energy of the solution 

keF r=0 

as computed from the full system (with M = 7 Legendre polynomials), the 
MZ reduced model with A = 2 without memory (keeping only the Markovian 
term) and the MZ reduced model with A = 2 with memory. Figure [2] shows 
the evolution of the standard deviation of the energy of the solution. The 
variance of the energy is given by 

I 1 

Var[E{t)] = ^ Yl Yl (27r)^'"fciri<ir-2^fe2r3<2r4t^nr-2r3r4-{E[-E(t)]}2, 
ki,k2&F T-i,...,r4=0 

where 

/■I 1 

drir2r3r4 = J ^ri (C)-^r2 (O-^rs (0-^r4 (C) 2 

The reduced model performs equally well with or without memory. Of 
course, the reduced model with memory is slower than the reduced model 
without memory. However, the reduced model with memory is still about 4 
times faster than the the full system. 

Figure [3] shows the evolution of the mean squared I2 norm of the gradient 
of the solution 

nG{t)] = Y,i2'^^k'\u,r{t)\'^^ 

k£F r=0 

as computed from the full system (with M = 7 Legendre polynomials), the 
MZ reduced model with A = 1 without memory (keeping only the Markovian 
term) and the MZ reduced model with A = 1 with memory. Figure S] shows 
the evolution of the standard deviation. The variance is given by 

1 

Var[G{t)] = E E (2^)^^1^2^fcin4ir2^fc2r3<2r4^nr2r3r4-{E[G(0]}^- 

ki,k2&F ri,...,r4=0 

It is obvious from the figures that the inclusion of the memory term im- 
proves the performance of the reduced model. Recall that the solution of 
Burgers equation is a contraction [lOj. Eventually, the complete description 
of the uncertainty caused by the uncertainty in the initial condition requires 
only a few polynomial chaos expansion coefficients. This happens at a time 
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Figure 3: Evolution of the mean of the squared I2 norm of the gradient of 
the ' ■■ ^ ^ ■ ^ ■ 




Figure 4: Evolution of the standard deviation of the squared I2 norm of 
the gradient of the solution calculated using only the first two Legendre 
polynomials. 
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scale that is dictated by the magnitude of the viscosity coefficient. That is 
why for long times the reduced model with and without memory have com- 
parable behavior to that of the full system. However, for short times, the 
inclusion of the memory term does make a difference because information 
from the higher chaos expansion coefficients is needed. The higher chaos 
expansion coefficients will have a more prolonged contribution for systems 
that possess unstable modes. In such cases, the inclusion of the memory 
term becomes imperative for short as well long times. Results for such cases 
will be presented elsewhere. 

5 Discussion and future work 

We have presented the application of the Mori-Zwanzig formalism to the 
construction of reduced models for systems of differential equations result- 
ing from polynomial chaos expansions of solutions of differential equations 
with initial condition uncertainty. In particular, we presented a way that 
the reduced model can be reformulated so that instead of integro-differential 
equations one has to solve differential equations. The problem that arises in 
any reduced model with memory is to compute the length of the memory. 
We have presented an algorithm which allows the estimation of the memory 
length on the fly. This means that one starts evolving the full system and 
use it to estimate the reduced model parameters. Once this is achieved, the 
simulation continues by evolving only the reduced model with the necessary 
parameters set equal to their estimated values from the first part of the algo- 
rithm. The construction presented here can also be applied to the problem 
of constructing reduced models for systems forced by random noise [9]. 

In the current work we have focus on the simplest form of the Markovian 
reformulation of the Mori-Zwanzig formalism which keeps only the first term 
in the expansion of the memory. In |14j , where we dealt with the problem of 
constructing reduced models when the original system exhibits parametric 
uncertainty we presented reduced models which used the first two terms in 
the memory expansion. However, the estimation of the necessary parame- 
ters for the reduced model was addressed in a different way. The method 
presented in |14j cannot be applied in the current case and this is the reason 
we had to devise an alternative approach. Results from the application of 
the current approach to reduced models where one keeps more terms in the 
memory expansion will be presented in a future publication. 

Finally, we mention that one can construct models which effect reduction 
both for the variables needed to describe uncertainty and the number of 
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variables needed to describe the system for one realization of the uncertainty 
sources. This two-level reduction is imperative in situations where solving 
even for one realization of the uncertainty sources is very expensive e.g. 
atmospheric flows, fluid structure interactions. 
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